/*In this implementation we use a multilayered perceptron for group feature selection */
  
#include<stdio.h>
#include<math.h>
#include<stdlib.h>

#define THRES 0.00001

float **data;
float **target;
int   *nFeaturesPerSensor;

float *output1;
float *output2;
float *input2;
float *output3;
float **weight12;
float **weight23;
float *beta;
float *previousBeta;
int *nFeaturesPerSensor;
int **featureSensorConnectivity;
int *sensorSignature;
   float *error;
   float *delta2;
   float *delta1;

int nDim;
int nData;
int nClass;

int final=0;

int nHiddenUnits;
int nSensors;
int Seed;

float eta=0.01;
float mu;

int iter;
int flagAtten=1;
int attenuationCount =0;

char dataFileName[40];

FILE *fpLogAtten;

void readData()
{ int i,j;
  FILE *fp;

   printf(" Enter the number of data points\n");
   scanf("%d",&nData);

   printf(" Enter the number of dimensions\n");
   scanf("%d",&nDim);

   printf(" Enter the number of classes\n");
   scanf("%d",&nClass);

   printf(" Enter the seed\n");
   scanf("%d",&Seed);


   printf(" Enter number of iterations \n");
   scanf("%d",&iter);

   printf("Enter the data file name:\n");
   scanf("%s",dataFileName);

   fp = fopen(dataFileName,"r");

     if(fp==NULL)

	 {  printf("cannot open file %s",dataFileName);

	    exit(0);
	 }

     data = (float**)malloc(nData*sizeof(float*));

	for( i = 0; i < nData; i++)
		{ data[i] = (float*)malloc(nDim*sizeof(float));

		   if(data[i]==NULL){ printf("Memory allocation failed\n");
				      exit(0);
				    }
		}
	target = (float**)malloc(nData*sizeof(float*));

	for( i = 0; i < nData; i++)
		{ target[i] = (float*)malloc(nClass*sizeof(float));

		   if(target[i]==NULL){ printf("Memory allocation failed\n");
					exit(0);
				      }
		}


   for(i=0;i<nData;i++)
     { for(j = 0;  j < nDim; j++) 
	  fscanf(fp,"%f",&data[i][j]);

       for(j = 0; j < nClass; j++)
	  fscanf(fp,"%f",&target[i][j]);
     }

   fclose(fp);

   printf("Data Read \n");

  printf("No of hidden units\n");
  scanf("%d",&nHiddenUnits);

  printf("No of sensors\n");
  scanf("%d",&nSensors);

  nFeaturesPerSensor = (int *) malloc(nSensors*sizeof(int));

/*   printf("enter number of features per sensor\n");
*/   
for(i =0; i<nSensors;i++)
nFeaturesPerSensor[i]=1;
 
/*    scanf("%d",&nFeaturesPerSensor[i]); */

  
 printf("enter eta\n");
 scanf("%f",&eta);
 printf("enter mu\n");
 scanf("%f",&mu);
}
 

void initialization()
{  int i,k,t,j;
 
  /* initializing weight matrix */

   weight23 = (float **)malloc(nClass*sizeof(float*));
   weight12 = (float **)malloc(nDim*sizeof(float*));

	for( i = 0; i < nClass; i++)
		{ weight23[i] = (float*)malloc(nHiddenUnits*sizeof(float));

		   if(weight23[i]==NULL){ printf("Memory allocation failed for weight\n");
					  exit(0);
					}
		}

	for( i = 0; i < nDim; i++)
		{ weight12[i] = (float*)malloc(nHiddenUnits*sizeof(float));

		   if(weight12[i]==NULL){ printf("Memory allocation failed for weight\n");
					  exit(0);
					}
		}

   srand((int)Seed);

   for(i = 0; i < nClass ; i++)
     for(k = 0; k < nHiddenUnits; k++)
	 weight23[i][k] = -0.5 + (float)drand48(); /*((rand()%987791)/987791.0) ;*/

   for(i = 0; i < nDim ; i++)
     for(k = 0; k < nHiddenUnits; k++)
	 weight12[i][k] = -0.5 + (float)drand48(); /* ((rand()%987791)/987791.0) ;*/

   /* allocating memory for output2 */
 
   output1  = (float *)malloc(nDim*sizeof(float*));

   output2  = (float *)malloc(nHiddenUnits*sizeof(float*));

   input2  = (float *)malloc(nHiddenUnits*sizeof(float*));

  /* allocating memory for output3 */

   output3 = (float *) malloc(nClass * sizeof(float));

  error = (float *)malloc(nClass * sizeof(float));



/* initializing feature connectivity */
 featureSensorConnectivity = (int **) malloc(nDim*sizeof(int *));

   for(i=0;i<nDim;i++)
     featureSensorConnectivity[i] = (int *) malloc(nSensors*sizeof(int));

  for(i=0;i<nDim;i++)
    for(j=0;j<nSensors;j++)
	 featureSensorConnectivity[i][j] = 0;

   sensorSignature = (int *)malloc(nDim*sizeof(int));

   for(i=0;i<nSensors;i++)
    { 
/** printf("the features related to sensor %d\n",i+1); **/

       for(j=0;j<nFeaturesPerSensor[i];j++)
	 { 
    /*** scanf("%d",&t); **/
	     t=i;
	     featureSensorConnectivity[t][i] = 1;
	     sensorSignature[t] = i;
	 }
    }

    

 /* for(i=0;i<nDim;i++)
    { for(j=0;j<nSensors;j++)
	printf("%d ",featureSensorConnectivity[i][j]);
	pirintf(" ");
	printf("\n");
    } */
  
    
	    

  /* initializing beta */

  beta = (float *) malloc(nSensors * sizeof(float));
  previousBeta = (float *) malloc(nSensors * sizeof(float));

   for(i = 0 ; i < nSensors; i++)
       { beta[i] = 2.5;
	 previousBeta[i]= 15;
       }

   delta2 = (float *) malloc(nHiddenUnits*sizeof(float));
   delta1 = (float *) malloc(nDim*sizeof(float));

 /* end initialization */

printf("initializationDone \n");
}

void resumePreviousSession()
{ int i,j,k,t;
  char sessionLogFile[40];
  FILE *fpLog,*fp;

  printf("Input the previous session log file on which you want to continue:\n");
  scanf("%s",sessionLogFile);

  fpLog = fopen(sessionLogFile,"r");
   if(fpLog==NULL)
     { printf("cannot open file %s \n",sessionLogFile);
       exit(0);
     }

  /* reading essential data */

   fscanf(fpLog,"%d",&nData);
   fscanf(fpLog,"%d",&nDim);
   fscanf(fpLog,"%d",&nClass);

   fscanf(fpLog,"%d",&iter);
   fscanf(fpLog,"%s",dataFileName);

   fp = fopen(dataFileName,"r");

     if(fp==NULL)

	 {  printf("cannot open file %s",dataFileName);
            fclose(fpLog);

	    exit(0);
	 }

     data = (float**)malloc(nData*sizeof(float*));

	for( i = 0; i < nData; i++)
		{ data[i] = (float*)malloc(nDim*sizeof(float));

		   if(data[i]==NULL){ printf("Memory allocation failed\n");
				      exit(0);
				    }
		}
	target = (float**)malloc(nData*sizeof(float*));

	for( i = 0; i < nData; i++)
		{ target[i] = (float*)malloc(nClass*sizeof(float));

		   if(target[i]==NULL){ printf("Memory allocation failed\n");
					exit(0);
				      }
		}


   for(i=0;i<nData;i++)
     { for(j = 0;  j < nDim; j++) 
	  fscanf(fp,"%f",&data[i][j]);

       for(j = 0; j < nClass; j++)
	  fscanf(fp,"%f",&target[i][j]);
     }

   fclose(fp);

   printf("Data Read \n");

  fscanf(fpLog,"%d",&nHiddenUnits);

  fscanf(fpLog,"%d",&nSensors);

  nFeaturesPerSensor = (int *) malloc(nSensors*sizeof(int));

/*   printf("enter number of features per sensor\n");
*/   
for(i =0; i<nSensors;i++)
nFeaturesPerSensor[i]=1;
 
/*    scanf("%d",&nFeaturesPerSensor[i]); */

  
 fscanf(fpLog,"%f",&eta);
 fscanf(fpLog,"%f",&mu);

 /* initializing weights */

   weight23 = (float **)malloc(nClass*sizeof(float*));
   weight12 = (float **)malloc(nDim*sizeof(float*));

	for( i = 0; i < nClass; i++)
		{ weight23[i] = (float*)malloc(nHiddenUnits*sizeof(float));

		   if(weight23[i]==NULL){ printf("Memory allocation failed for weight\n");
					  exit(0);
					}
		}

	for( i = 0; i < nDim; i++)
		{ weight12[i] = (float*)malloc(nHiddenUnits*sizeof(float));

		   if(weight12[i]==NULL){ printf("Memory allocation failed for weight\n");
					  exit(0);
					}
		}


   for(i = 0; i < nClass ; i++)
     for(k = 0; k < nHiddenUnits; k++)
	 fscanf(fpLog,"%f",&weight23[i][k]) ;

   for(i = 0; i < nDim ; i++)
     for(k = 0; k < nHiddenUnits; k++)
	 fscanf(fpLog,"%f",&weight12[i][k]);

   /* allocating memory for output2 */
 
   output1  = (float *)malloc(nDim*sizeof(float*));

   output2  = (float *)malloc(nHiddenUnits*sizeof(float*));

   input2  = (float *)malloc(nHiddenUnits*sizeof(float*));

  /* allocating memory for output3 */

   output3 = (float *) malloc(nClass * sizeof(float));

  error = (float *)malloc(nClass * sizeof(float));



/* initializing feature connectivity */
 featureSensorConnectivity = (int **) malloc(nDim*sizeof(int *));

   for(i=0;i<nDim;i++)
     featureSensorConnectivity[i] = (int *) malloc(nSensors*sizeof(int));

  for(i=0;i<nDim;i++)
    for(j=0;j<nSensors;j++)
	 featureSensorConnectivity[i][j] = 0;

   sensorSignature = (int *)malloc(nDim*sizeof(int));

   for(i=0;i<nSensors;i++)
    { 
/** printf("the features related to sensor %d\n",i+1); **/

       for(j=0;j<nFeaturesPerSensor[i];j++)
	 { 
    /*** scanf("%d",&t); **/
	     t=i;
	     featureSensorConnectivity[t][i] = 1;
	     sensorSignature[t] = i;
	 }
    }

    


  /* initializing beta */

  beta = (float *) malloc(nSensors * sizeof(float));
  previousBeta = (float *) malloc(nSensors * sizeof(float));

   for(i = 0 ; i < nSensors; i++)
        fscanf(fpLog,"%f",&beta[i]);
       

   delta2 = (float *) malloc(nHiddenUnits*sizeof(float));
   delta1 = (float *) malloc(nDim*sizeof(float));

 /* end initialization */

printf("initializationDone \n");

fclose(fpLog);

}
 

void saveCurrentSession()
{ int i,j,k,t;
  char sessionLogFile[40];
  FILE *fpLog,*fp;

  printf("Input the session log file on which you want to save:\n");
  scanf("%s",sessionLogFile);

  fpLog = fopen(sessionLogFile,"w");
   if(fpLog==NULL)
     { printf("cannot open file %s \n",sessionLogFile);
       exit(0);
     }

  /* writing  essential data */

   fprintf(fpLog,"%d\n",nData);
   fprintf(fpLog,"%d\n",nDim);
   fprintf(fpLog,"%d\n",nClass);

   fprintf(fpLog,"%d\n",iter);
   fprintf(fpLog,"%s\n",dataFileName);

   fprintf(fpLog,"%d\n",nHiddenUnits);
   fprintf(fpLog,"%d\n",nSensors);

  
 fprintf(fpLog,"%f\n",eta);
 fprintf(fpLog,"%f\n",mu);

 /* writing weights */



   for(i = 0; i < nClass ; i++)
     { for(k = 0; k < nHiddenUnits; k++)
	 fprintf(fpLog,"%f ",weight23[i][k]) ;
         fprintf(fpLog,"\n");
     }

   for(i = 0; i < nDim ; i++)
     { for(k = 0; k < nHiddenUnits; k++)
	 fprintf(fpLog,"%f ",weight12[i][k]);
         fprintf(fpLog,"\n");
     }

  /* writing beta */

   for(i = 0 ; i < nSensors; i++)
        fprintf(fpLog,"%f ",beta[i]);
         fprintf(fpLog,"\n");

printf("log written\n");

fclose(fpLog);

}
 



float calculateOutput1(dimID, dataID)
int dimID;
int dataID;
{ int i,j;
  float temp,temp1;

   temp = (float)pow(beta[sensorSignature[dimID]], 2.0);
   
   temp1 =(float) exp(-temp);

    temp = temp1 * data[dataID][dimID];

    /* if(temp==0) temp =0.0001;*/
    return(temp);
}


float calculateOutput2( hiddenUnitID, dataID)
 int hiddenUnitID;
 int dataID;

{ int i,j;
  float out,temp =0;

   for(i=0;i<nDim;i++)
      temp = temp + weight12[i][hiddenUnitID]* output1[i];

   out = 1/(1 + (float)exp(-temp));

  return(out);
}



float calculateOutput3(classID)
int classID;
{ int i,k;
  float out=0,out1;


    for( k = 0; k < nHiddenUnits ; k++)
      out = out + weight23[classID][k] * output2[k];

  out1 = 1/(1 + (float)exp(-out));

  return(out1);
} 


forward(dataNo)
int dataNo;

{ int i,j,l,k;
  float temp;

      for(k=0; k< nDim ; k++)
	   output1[k] = calculateOutput1(k,dataNo);

      for (k = 0; k < nHiddenUnits; k++)
	   output2[k] = calculateOutput2(k,dataNo);
      

	for( i = 0; i < nClass; i++)
	     output3[i] = calculateOutput3(i);
}

checkAttenuationTermination()
{ int i,j;
  float temp=0;

   for(i=0;i<nSensors;i++)
     temp += pow(beta[i] - previousBeta[i],2.0);

    temp = (float)sqrt((double)temp);

    if(temp<THRES) attenuationCount++;
	else attenuationCount =0;

     if(attenuationCount==500) { flagAtten =0;
				 printf("attenuation terminated \n");
			       }

   for(i=0;i<nSensors;i++)
	previousBeta[i] = beta[i];
    

}

void learning()
{  int i,j,k,l,mis;
   int iterNo=0,dataID;
   float squareError;
   float extra;


 while(iterNo < iter)
  { for( dataID = 0; dataID < nData; dataID++)
     { forward(dataID);

       for( k = 0; k < nHiddenUnits; k++)
	    delta2[k] = 0.0;

       for( k = 0; k < nDim; k++)
	    delta1[k] = 0.0;

	for( i = 0; i < nClass; i++)
	    error[i] = output3[i] - target[dataID][i];
           
    
       /* Weight updation */

	for( i = 0; i < nClass; i++)
	    for( k = 0; k < nHiddenUnits; k++)
		weight23[i][k] = weight23[i][k] - eta * error[i] * output2[k]*output3[i]*(1-output3[i]);


	   for(k = 0; k < nHiddenUnits; k++)
		for( i = 0; i < nClass; i++) 
		      delta2[k] +=  error[i] * weight23[i][k]*output3[i]*(1-output3[i]);


	    for(i=0; i< nDim; i++)
	       for(j=0;j<nHiddenUnits;j++)
		 weight12[i][j] = weight12[i][j] - eta * delta2[j]*output1[i]*(1-output2[j])*output2[j];

		 /*  attenuation updation */ 
 if(flagAtten==1)
  {
	  for(k=0;k<nDim;k++)
	     for(i=0;i<nHiddenUnits;i++)
		 delta1[k] += delta2[i] * weight12[k][i] * (1-output2[i]) * output2[i];

	  extra=0;
 
	      for( i=0; i<nSensors; i++)
		  { for(k=0; k< nDim; k++)
		       { if(featureSensorConnectivity[k][i]==1)
			   extra = extra - 2*beta[i]*delta1[k]*output1[k];
		       }
		     beta[i] = beta[i] - mu * extra; 
		     extra =0;
		  /*     beta[i] = beta[i] - mu * (extra - (nClass/(float)nSensors)*beta[i]*exp(-beta[i]*beta[i])); */
		  }
   }
     
       
     }

  squareError =0;

 /* if(flagAtten==1) checkAttenuationTermination();*/
     
      mis = misClass();
    for(i=0; i<nSensors; i++)  
     fprintf(fpLogAtten,"%f ",(float)exp(-beta[i]*beta[i]));
     fprintf(fpLogAtten,"%d\n",mis);

  if(iterNo%50==0)
   { 
    for( dataID = 0; dataID < nData; dataID++)
     { forward(dataID);
	 for(i = 0; i < nClass; i++)
	    squareError = squareError + (output3[i]-target[dataID][i]) *(output3[i]-target[dataID][i]);
     }
       printf("SSE after iteration %d is %f\n",iterNo,squareError);
       squareError =0;

      mis = misClass();
     if(mis== 0)return;

  printf("Misclassification = %d\n", mis);
 
if ((iter-1)==iterNo)    for(i=0; i<nSensors; i++)  
     	printf("%f ",(float)exp(-beta[i]*beta[i]));
	/*changed to above-22Dec04*printf("%f\n",(float)exp(-beta[i]*beta[i]));*/
    
/*for(i=0; i<nSensors; i++)  
     fprintf(fpLogAtten,"%f ",(float)exp(-beta[i]*beta[i]));
     fprintf(fpLogAtten,"%d\n",mis); */

     /* printf("%f\n",beta[i]); */

    }

   iterNo++;

  }
}

 int misClass()    
{ int i,j,k, nMisClass=0, label=0,dataNo;
  float max=0;
  
  /*** changes made by NRP for confusion table ***/
  int confusion[30][30];
  if (final==1)
  { for (i=0;i<30;i++)
    for (j=0;j<30;j++)
       confusion[i][j]=0;
   }
/*** Changes end here *****/

   for(dataNo=0; dataNo < nData; dataNo++)
     { 
	forward(dataNo);
	for(i=0; i< nClass; i++)
	  {
		if(output3[i]>max)
		   { 
			max = output3[i];
			      label = i;
		   }
	  }
 
      if(target[dataNo][label]!=1) nMisClass++;
      max =0;

      
     
      /**** Changes made by NRP ****/
      
    if (final==0) label =0;
    
    if(final==1)    
      {   j=0;
	  for (i=0; i< nClass; i++)
	    if (target[dataNo][i]==1) j=i;
	    
	     confusion[j][label]++;
	     label=0;
       }
      /****** Changes finish here **/
   }

 /* printf("Misclassification = %d\n", nMisClass);*/
 if (final==1)
   {
	for (i=0;i<nClass;i++)
	 {
		for (j=0;j<nClass;j++)
		  printf("%d ",confusion[i][j]);
		printf("\n");
	 }
   }
  return(nMisClass);
}

void toTest()
{ int i,j,k,mis;
  char choice[10], fileNam[30];
  FILE *fpTest;

  printf("Want to test?yes/no\n");
  scanf("%s",choice);

  if(strcmp(choice,"yes")==0)
  {  printf("input testDataFile\n");
     scanf("%s",fileNam);
     
     fpTest = fopen(fileNam,"r");

      if(fpTest==NULL){printf("cannot open file %s\n",fileNam);
		       exit(0);
		      }
     for(i=0;i<nData;i++)
      { free(data[i]);
	free(target[i]);
      }
       free(data);
       free(target);
    
    printf("enter no of test data points\n");
    scanf("%d",&nData);
    printf("scanned\n");

    printf("Allocating memory----");
	   
       data = (float**)malloc(nData*sizeof(float*));

	for( i = 0; i < nData; i++)
		{ data[i] = (float*)malloc(nDim*sizeof(float));

		   if(data[i]==NULL){ printf("Memory allocation failed\n");
				      exit(0);
				    }
		}
	target = (float**)malloc(nData*sizeof(float*));

	for( i = 0; i < nData; i++)
		{ target[i] = (float*)malloc(nClass*sizeof(float));

		   if(target[i]==NULL){ printf("Memory allocation failed\n");
					exit(0);
				      }
		}
     printf("allocated\n");

printf("reading data ---- ");
   for(i=0;i<nData;i++)
     { for(j = 0;  j < nDim; j++)
	  fscanf(fpTest,"%f",&data[i][j]);

       for(j = 0; j < nClass; j++)
	  fscanf(fpTest,"%f",&target[i][j]);
     }
   fclose(fpTest);
   printf("read\n");
   mis =misClass();
   printf("Misclassification on test data = %d",mis);
  }
}
    
main()
{ int i,j,k,choice;
  char choiceString[10];
  FILE *fp;

  printf("Options:\n");
  printf("1) resume a previous session\n");
  printf("2) start a fresh session\n\n");
  printf("enter Choice: ");
  scanf("%d",&choice);
  printf("\n");

  if(choice == 1)
    resumePreviousSession();
  else if( choice == 2)
     { readData();
       initialization();
     }
  else  { printf("invalid choice -- session terminated, start again\n");
          exit(0);
        }


  fpLogAtten = fopen("attenuationLog1.log","w");


  learning();
  fclose(fpLogAtten);
  printf("training terminated\n");
  final=1;
  printf("Misclassification on training Data = %d\n",misClass());

  printf("Want to save current session?yes/no\n");
  scanf("%s",choiceString);

  if(strcmp(choiceString,"yes")==0)
     saveCurrentSession(); 
  else printf("session not saved\n");

fp = fopen("finalAttenuators.dat","w");

  for(i=0; i<nSensors; i++)
     fprintf(fp,"%f\n",(float)exp(-beta[i]*beta[i]));
 fclose(fp);
  
   toTest();
}  
